home *** CD-ROM | disk | FTP | other *** search
/ Compendium Deluxe 2 / LSD and 17bit Compendium Deluxe - Volume II.iso / a / prog / cprog / illum.lha / D.c < prev    next >
Encoding:
C/C++ Source or Header  |  1988-05-11  |  14.0 KB  |  420 lines
  1. /* ****************************************************************
  2.  *                            D.c
  3.  * ****************************************************************
  4.  *  MODULE PURPOSE:
  5.  *    This module contains microfacet distribution routines and
  6.  *    various support routines.  One of support functions includes
  7.  *    computing the coefficients corresponding to beta for each
  8.  *    function (as described in Section 4.3.1-Improved Specular
  9.  *    Reflection).  Beta is the angle between H and N where the
  10.  *    function is equal to half the value at N = H.  Beta is in
  11.  *    radians.
  12.  *
  13.  *  MODULE CONTENTS:
  14.  *    D_phong_init    - compute Ns given beta
  15.  *    D_phong         - Phong cosine power function
  16.  *    D_blinn_init    - compute Ns given beta
  17.  *    D_blinn         - Blinn cosine power function
  18.  *    D_gaussian_init - compute C1 given beta
  19.  *    D_gaussian      - Gaussian distribution
  20.  *    D_reitz_init    - compute c2 given beta
  21.  *    D_reitz         - Trowbridge-Reitz
  22.  *    D_cook_init     - compute m given beta
  23.  *    D_cook          - Cook model
  24.  *    D_g             - compute apparent roughness
  25.  *    D_tau           - compute correlation distance
  26.  *    D_sigma         - compute rms roughness
  27.  *    D_m             - compute slope
  28.  *    D_coherent      - coherent roughness attenuation
  29.  *                        (per Beckmann)
  30.  *    D_incoherent    - incoherent microfacet attenuation
  31.  *                        (per Beckmann and Bahar)
  32.  *    D_Vxy           - Vxy used in D_incoherent()
  33.  *
  34.  *  ASSUMPTIONS:
  35.  *  > The following are defined in math.h
  36.  *          HUGE        largest floating pont value
  37.  *          M_PI        pi
  38.  *          M_PI_2      pi/2
  39.  *          M_PI_4      pi/4
  40.  *          M_SQRT2     root of 2
  41.  *
  42.  *  > The direction vectors N, V, and L are assumed to be oriented
  43.  *      to the same side of the surface.
  44.  */
  45. #include <stdio.h>
  46. #include <math.h>
  47. #include "geo.h"
  48. #include "D.h"
  49.  
  50. /* ****************************************************************
  51.  * double D_phong_init (beta)
  52.  *  double      beta        (in) - angle between N and
  53.  *                              H where function = .5
  54.  *  Returns Ns given beta for the Phong (1975) specular function,
  55.  *   per footnote to Eq.(\*(rH)
  56.  */
  57. double D_phong_init (beta)
  58. double          beta;
  59. {   if (beta <= 0.0) return HUGE;
  60.     if (beta >= M_PI_4) return 0.0;
  61.     return -(log(2.0) / log(cos(2.0 * beta)));
  62. }
  63.  
  64. /* ****************************************************************
  65.  * double D_phong (N, L, V, Ns)
  66.  *  DIR_VECT    *N, *L, *V  (in) - N, L, V vectors
  67.  *  double      Ns          (in) - specular exponent
  68.  *
  69.  *  Returns the value of the Phong (1975) specular function given
  70.  *   surface normal, incident light direction, view direction, and
  71.  *   the specular exponent.  Refer to Eqs.(\*(rU) and (\*(rI).
  72.  */
  73. double D_phong (N, L, V, Ns)
  74. DIR_VECT        *N, *L, *V;
  75. double          Ns;
  76. {   double      Rv_dot_L;
  77.     if ((Rv_dot_L = geo_dot(geo_rfl(V,N),L)) < 0.0) return 0.0;
  78.     return pow (Rv_dot_L, Ns);
  79. }
  80.  
  81. /* ****************************************************************
  82.  * double D_blinn_init (beta)
  83.  *  double      beta        (in) - angle between N and
  84.  *                              H where function = .5
  85.  *  Returns Ns given beta for the cosine power distribution
  86.  *   function presented by Blinn (1977).  Refer to Eq.(\*(rHa)
  87.  */
  88. double D_blinn_init (beta)
  89. double          beta;
  90. {   if (beta <= 0.0) return HUGE;
  91.     if (beta >= M_PI_2) return 0.0;
  92.     return -(log(2.0) / log(cos(beta)));
  93. }
  94.  
  95. /* ****************************************************************
  96.  * double D_blinn (N, L, V, Ns)
  97.  *  DIR_VECT    *N, *L, *V  (in) - N, L, V vectors
  98.  *  double      Ns          (in) - specular exponent
  99.  *
  100.  *  Returns the cosine power distribution function presented by
  101.  *   Blinn (1977) given surface normal, incident light direction,
  102.  *   view direction, and the specular exponent. Refer to Eq.(\*(e5)
  103.  */
  104. double D_blinn (N, L, V, Ns)
  105. DIR_VECT        *N, *L, *V;
  106. double          Ns;
  107. {   double      N_dot_H;
  108.     if ((N_dot_H = geo_dot(N,geo_H(V,L))) < 0.0) return 0.0;
  109.     return pow (N_dot_H, Ns);
  110. }
  111.  
  112. /* ****************************************************************
  113.  * double D_gaussian_init (beta)
  114.  *  double      beta        (in) - angle between N and
  115.  *                              H where function = .5
  116.  *  Returns C1 given beta for the Gaussian distribution presented
  117.  *   by Blinn (1977).  Refer to Eq.(\*(rHb)
  118.  */
  119. double D_gaussian_init (beta)
  120. double          beta;
  121. {   if (beta <= 0.0) return HUGE;
  122.     return sqrt(log(2.0)) / beta;
  123. }
  124.  
  125. /* ****************************************************************
  126.  * double D_gaussian (N, L, V, C1)
  127.  *  DIR_VECT    *N, *L, *V  (in) - N, L, V vectors
  128.  *  double      C1          (in) - shape coefficient
  129.  *
  130.  *  Returns the Gaussian distribution function presented by
  131.  *   Blinn (1977) given surface normal, incident light direction,
  132.  *   view direction, and the specular exponent. Refer to Eq.(\*(eA)
  133.  */
  134. double D_gaussian (N, L, V, C1)
  135. DIR_VECT        *N, *L, *V;
  136. double          C1;
  137. {   double      tmp;
  138.     double      N_dot_H;
  139.     if ((N_dot_H = geo_dot(N,geo_H(V,L))) < 0.0) return 0.0;
  140.     tmp = acos(N_dot_H) * C1;
  141.     return exp (-(tmp * tmp));
  142. }
  143.  
  144. /* ****************************************************************
  145.  * double D_reitz_init (beta)
  146.  *  double      beta        (in) - angle between N and
  147.  *                              H where function = .5
  148.  *  Returns C2 squared given beta for the Trowbridge-Reitz
  149.  *   distribution function presented by Blinn (1977). Refer to
  150.  *   Eq.(\*(rHc). C2 is squared for computational efficiency later.
  151.  */
  152. double D_reitz_init (beta)
  153. double          beta;
  154. {   double      cos_beta;
  155.     if (beta <= 0.0) return 0.0;
  156.     cos_beta = cos(beta);
  157.     return ((cos_beta * cos_beta) - 1.0) /
  158.     ((cos_beta * cos_beta) - sqrt(2.0));
  159. }
  160.  
  161. /* ****************************************************************
  162.  * double D_reitz (N, L, V, C2_2)
  163.  *  DIR_VECT    *N, *L, *V  (in) - N, L, V vectors
  164.  *  double      C2_2        (in) - C2 squared
  165.  *
  166.  *  Returns the Trowbridge-Reitz distribution function presented
  167.  *   by Blinn (1977) given surface normal, incident light
  168.  *   direction, view direction, and the specular exponent. Refer
  169.  *   to Eq.(\*(eB)
  170.  */
  171. double D_reitz (N, L, V, C2_2)
  172. DIR_VECT        *N, *L, *V;
  173. double          C2_2;
  174. {   double      tmp;
  175.     double      N_dot_H;
  176.     if ((N_dot_H = geo_dot(N,geo_H(V,L))) < 0.0) return 0.0;
  177.     tmp = C2_2 / ((N_dot_H * N_dot_H * (C2_2 - 1.0)) + 1.0);
  178.     return tmp * tmp;
  179. }
  180.  
  181. /* ****************************************************************
  182.  * double D_cook_init (beta)
  183.  *  double      beta        (in) - angle between N and
  184.  *                              H where function = .5
  185.  *  Returns m squared corresponding to beta, refer to Eq.(\*(rT).
  186.  *   m is squared for computational efficiency in later use.
  187.  */
  188. double D_cook_init (beta)
  189. double          beta;
  190. {   double      tan_beta;
  191.     if (beta <= 0.0) return 0.0;
  192.     if (beta >= M_PI_2) return HUGE;
  193.     tan_beta = tan(beta);
  194.     return -(tan_beta * tan_beta) /
  195.     log(pow(cos(beta),4.0)/2.0);
  196. }
  197.  
  198. /* ****************************************************************
  199.  * double D_cook (N, L, V, m_2)
  200.  *  DIR_VECT    *N, *L, *V  (in) - N, L, V vectors
  201.  *  double      m_2         (in) - m squared
  202.  *
  203.  *  Returns microfacet distribution probability (Cook 1982) derived
  204.  *   from (Beckmann 1963).  Refer to Eq.(\*(rX).
  205.  */
  206. double D_cook (N, L, V, m_2)
  207. DIR_VECT        *N, *L, *V;
  208. double          m_2;
  209. {   double      tmp;
  210.     double      N_dot_H;
  211.     if ((N_dot_H = geo_dot(N,geo_H(V,L))) < 0.0) return 0.0;
  212.     tmp = -(1.0 - (N_dot_H * N_dot_H)) / (N_dot_H * N_dot_H * m_2);
  213.     return exp(tmp)/(4.0 * M_PI * m_2 * pow(N_dot_H,4.0));
  214. }
  215.  
  216. /* ****************************************************************
  217.  * double D_g (N, L, V, sigma, lambda)
  218.  *  DIR_VECT    *N, *L, *V  (in) - N, L, V vectors
  219.  *  double      sigma       (in) RMS roughness, nm
  220.  *  double      lambda      (in) wavelength, nm
  221.  *
  222.  *  Returns the apparent roughness, g, as given by Beckmann (1963).
  223.  *   Refer to Eq.(\*(rY).
  224.  */
  225. double D_g (N, L, V, sigma, lambda)
  226. DIR_VECT        *N, *L, *V;
  227. double          sigma, lambda;
  228. {   double      tmp;
  229.     tmp = geo_dot(N,L) + geo_dot(N,V);
  230.     return (4.0 * M_PI * M_PI * sigma * sigma * tmp * tmp) /
  231.     (lambda * lambda);
  232. }
  233.  
  234. /* ****************************************************************
  235.  * double D_tau (sigma, m)
  236.  *  double      sigma       (in) rms roughness (nm)
  237.  *  double      m           (in) rms slope
  238.  *
  239.  *  Returns tau (correlation distance) in nm.  The relationship of
  240.  *   roughness, slope, and correlation distance given by
  241.  *   Beckmann (1963) is assumed, refer to Eq.(\*(rD).
  242.  */
  243. double D_tau (sigma, m)
  244. double          sigma, m;
  245. {   return (2.0 * sigma) / m;
  246. }
  247.  
  248. /* ****************************************************************
  249.  * double D_sigma (m, tau)
  250.  *  double      m           (in) rms slope
  251.  *  double      tau         (in) correlation distance (nm)
  252.  *
  253.  *  Returns sigma (rms roughness). The relationship of roughness,
  254.  *   slope, and correlation distance given by Beckmann (1963) is
  255.  *   assumed, refer to Eq.(\*(rD).
  256.  */
  257. double D_sigma (m, tau)
  258. double          m, tau;
  259. {   return (m * tau)/ 2.0;
  260. }
  261.  
  262. /* ****************************************************************
  263.  * double D_m (sigma, tau)
  264.  *  double      sigma       (in) rms roughness (nm)
  265.  *  double      tau         (in) correlation distance (nm)
  266.  *
  267.  *  Returns m (rms slope). The relationship  of roughness, slope,
  268.  *   and correlation distance given by Beckmann (1963) is assumed,
  269.  *   refer to Eq.(\*(rD).
  270.  */
  271. double D_m (sigma, tau)
  272. double          sigma, tau;
  273. {   return (2.0 * sigma) / tau;
  274. }
  275.  
  276. /* ****************************************************************
  277.  * double D_coherent (g)
  278.  *  double      g           (in) apparent roughness
  279.  *
  280.  *  Returns the coherent roughness attenuation given
  281.  *  by Beckmann (1963), refer to Eq.(\*(rZ).
  282.  */
  283. double D_coherent (g)
  284. double g;
  285. {   return exp(-g);
  286. }
  287.  
  288. /* ****************************************************************
  289.  * double D_incoherent (N, L, V, m, g, lambda, tau)
  290.  *  DIR_VECT    *N, *L, *V  (in) N, L, V vectors
  291.  *  double      m           (in) m (slope)
  292.  *  double      g           (in) apparent roughness
  293.  *  double      lambda      (in) wavelength (nm)
  294.  *  double      tau         (in) correlation distance (nm)
  295.  *
  296.  *  Returns the microfacet distribution function given by
  297.  *   Beckmann (1963) as interpreted by Koestner (1986), refer to
  298.  *   Eq.(\*(rC).
  299.  */
  300. double D_incoherent (N, L, V, m, g, lambda, tau)
  301. DIR_VECT        *N, *L, *V;
  302. double          m, g, lambda, tau;
  303. {   DIR_VECT    *H;
  304.     double      D1, D2;
  305.     double      denom;
  306.     double      N_dot_H, N_dot_H_2;
  307.  
  308.     if (g <= 0.0) return 0.0;
  309.     if ((H = geo_H(V,L)) == NULL) return 0.0;
  310.     if ((N_dot_H = geo_dot(N,H)) <= 0.0) return 0.0;
  311.     N_dot_H_2 = N_dot_H * N_dot_H;
  312.     denom = 4.0 * M_PI * m * m * N_dot_H_2 * N_dot_H_2;
  313.  
  314.     /* if g < 8.0, evaluate the series expansion, Eq.(\*(rCb)
  315.      *  terminating when a term falls below 0.00001.  For small g
  316.      *  the evaluation terminates quickly.  For large g the series
  317.      *  converges slowly.  If g is near 8.0 the first terms of the
  318.      *  series will be less that the termination criteria.  This
  319.      *  requires an additional termination criteria that the values
  320.      *  of the terms are decreasing at the time the termination
  321.      *  criteria is reached.
  322.      */
  323.     if (g < 8.0) {
  324.     double      inc, ct, ct_factorial, g_pow;
  325.     double      last_inc, Vxz_t_over_4;
  326.     Vxz_t_over_4 = (tau * tau * D_Vxz(N,L,V,lambda)) / 4.0;
  327.     D1 = 0.0;
  328.     ct = 1.0;
  329.     ct_factorial = 1.0;
  330.     g_pow = g * g;
  331.     inc = 0.0;
  332.     do {
  333.         last_inc = inc;
  334.         inc = (g_pow / (ct * ct_factorial)) *
  335.         exp(-(g + Vxz_t_over_4/ct));
  336.         D1 += inc;
  337.         ct += 1.0;
  338.         ct_factorial *= ct;
  339.         g_pow *= g;
  340.     } while (((inc/denom) > 0.00001) || (inc > last_inc));
  341.     D1 /= denom;
  342.  
  343.     /* if g < 5.0, only the series expansion is used.  If
  344.      *  5.0 < g < 8.0, then we are in a transition range
  345.      *  for interpolation between the series expansion and
  346.      *  the convergence expression.
  347.      */
  348.     if (g < 5.0) return D1;
  349.     }
  350.  
  351.     /* if g > 5.0, evaluate the convergence expression, Eq.(\*(rCc)
  352.      *  (this routine would have returned earlier if G < 5.0).
  353.      */
  354.     {   double  tmp;
  355.     tmp = (N_dot_H_2 - 1.0) / (N_dot_H_2 * m * m);
  356.     D2 = exp(tmp) / denom;
  357.  
  358.     /* if g > 8.0, only the convergence expression is used,
  359.      *  otherwise we are in a transition zone between the
  360.      *  convergent expression and series expansion.
  361.      */
  362.     if (g > 8.0) return D2;
  363.     }
  364.  
  365.     /* linear interpolation between D1 and D2 for
  366.      *          5.0 < g < 8.0
  367.      */
  368.     return (((8.0 - g) * D1) + ((g - 5.0) * D2)) / 3.0;
  369. }
  370.  
  371. /* ****************************************************************
  372.  * double D_Vxz (N, L, V, lambda)
  373.  *  DIR_VECT    *N, *L, *V  (in) N, L, V vectors
  374.  *  double      lambda      (in) wavelength (nm)
  375.  *
  376.  *  Returns the value of Vxz per Eq.(\*(rCd).
  377.  */
  378. double D_Vxz (N, L, V, lambda)
  379. DIR_VECT        *N, *L, *V;
  380. double          lambda;
  381. {   double      N_dot_L, N_dot_V, cos_phi;
  382.     DIR_VECT    sin_i, sin_r;
  383.     double      len_2_i, len_2_r;
  384.  
  385.     /* compute the sine squared of the incident angle
  386.      */
  387.     N_dot_L = geo_dot(N, L);
  388.     len_2_i = 1.0 - (N_dot_L * N_dot_L);
  389.  
  390.     /* compute the sine squared of the reflected angle
  391.      */
  392.     N_dot_V = geo_dot(N, V);
  393.     len_2_r = 1.0 - (N_dot_V * N_dot_V);
  394.  
  395.     /* if the sine of either the incident or the
  396.      *  reflected angle is zero, then the middle
  397.      *  term is zero.
  398.      */
  399.     if ((len_2_i > 0.0) && (len_2_r > 0.0)) {
  400.     /* the two sine vectors are of length equal to the sine of
  401.      *  the angles.  The dot product is the cosine times the
  402.      *  lengths of the vectors (sines of the angles).
  403.      */
  404.     sin_i.i = L->i - (N_dot_L * N->i);
  405.     sin_i.j = L->j - (N_dot_L * N->j);
  406.     sin_i.k = L->k - (N_dot_L * N->k);
  407.     sin_r.i = V->i - (N_dot_V * N->i);
  408.     sin_r.j = V->j - (N_dot_V * N->j);
  409.     sin_r.k = V->k - (N_dot_V * N->k);
  410.     cos_phi = geo_dot (&sin_i, &sin_r);
  411.     return (4.0 * M_PI * M_PI * (len_2_i +
  412.         (2.0 * cos_phi) + len_2_r)) / (lambda * lambda);
  413.     }
  414.     else {
  415.     return (4.0 * M_PI * M_PI * (len_2_i + len_2_r)) /
  416.         (lambda * lambda);
  417.     }
  418. }
  419. /* ************************************************************* */
  420.